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ABSTRACT 


We have investigated the effect of stellar encounters on the formation and 
disruption of the Oort cloud using the classical impulse approximation. We cal¬ 
culate the evolution of a planetesimal disk into a spherical Oort cloud due to 
the perturbation from passing stars for 10 Gyr. We obtain the empirical fits 
of the e-folding time for the number of Oort cloud comets using the standard 
exponential and Kohlrausch formulae as functions of the stellar parameters and 
the initial semimajor axes of planetesimals. The e-folding time and the evolution 
timescales of the orbital elements are also analytically derived. In some calcula¬ 
tions, the effect of the Galactic tide is additionally considered. We also show the 
radial variations of the e-folding times to the Oort cloud. From these timescales, 


we show that if the initial planetesimal disk has the semim ajor axes distri 


Dution 


dn/da oc a which is produced by planetary scattering flHiguchi et ah 


200611 . 


the e-folding time for planetesimals in the Oort cloud is ~10 Gyr at any heliocen¬ 
tric distance r. This uniform e-folding time over the Oort cloud means that the 
supply of comets from the inner Oort cloud to the outer Oort cloud is sufficiently 
effective to keep the comet distribution as dn/dr oc r“^. We also show that 
the dual distribution of the semimajor axes in the Oort cloud is approximately 
proportional to for any initial distribution. 


Subject headings: Oort Cloud — comets: general 
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INTRODUCTION 


_ ^The Oort Cloud is a spherical comet reservoir surrounding the Solar system fjOort 

1950). Observations and statistical studies estimate t hat it consists of m ore than 10^^ 


comets and is on the order of 10^-10® AU in size (e.g., 


Pones et al 


20041) ■ This structure 


can be described as an assembly of comets whose perihelion distances extend outside the 
planetary region, aphelion distances smaller than the tidal radius of the Sun (~1 pc), and a 
nearly isotropic inclination distribution. Now it is generally accepted that the comets are 
residual planetesimals from planet formation and are originally inside the planetary region. 

In the standard scenario of the Oort cloud formation, these residual planetesimals with 
small eccentricities and inclinations are scattered by giant planets and their semimajor axes 
and eccentricities are raised. Such planetesimals with large aphelion distances are affected 
by external forces. The perihelion distances are pulled out of the planetary region and the 
inclinations are randomized by external forces and then the spherical structure of the Oort 
cloud is attained. The perturbation from passing stars is the only external force considered 
in the original Oort scenario. Many authors have studied the effect of the external forces 
and now we recognize that not only passing s tars, but also the Galactic disk and giant 


molecular clouds are effective perturbers (e.g.. 


Pones et al. 


2nn4 . 


Among these ex 

ernal 

Harrington 

19^5; 

By] 

1986: 


brces, the most effective one is the Galactic tide (e.g.. 


Heisler fc Tremaine 


19861 ). The Galactic tide can efficiently 


raise the perihelion distances of planetesimals with large semimajor axes and eccentricities. 


Higuchi et al. 


(120071) derived the analytical formulae for time evolution of orbital elements of 


planetesimals under the effect of the Galactic tide. They considered the vertical component 
of the Galactic tide against the Galactic plane, which is ~ 10 times larger than the other 
components. They showed that the inclination distribution attained by the Galactic tide in 
5 Gyr is far from the isotropic one, which is expected from the observations of long-period 




































4 


comets. This is because the Galactic tide causes, not the randomization, but the periodic 
oscillations of orbital elements and the synchronized evolution of the orbital elements 
produces two peaks in the inclination distribution. Also, the inclination range attained by 
the vertical component of the Galactic tide is not from 0 to 180°, but the function of the 
inclination of the ecliptic plane against the Galactic plane. In the case of the Solar system, 
a planetesimal initially on the ecliptic plane cannot attain an inclination larger than 153°. 


Passing stars randomize and eject planetesimals from the Solar system by giving 
velocity kicks, while the vertical component of the Galactic tidal force does not change the 
energy of the planetesimals. The perturbations from passing stars, due to their random-walk 
nature, may play an important role in the production of the nearly isotropic inclination 
distribution of planetesimals. Many authors have examined this effect, main ly by analytical 


approach (e.g.. 


Heisler et al. 


1987 


Dybczyhski 


2002 


Fouchard et al. 


2011). Their main 


interest is in the production of long-period comets from the spherical Oort cloud. 


The direct simulation of the Oort cloud formation including th e perturbation s froni 


giant planets, 
and redone by 


he Galactic c 


Pones et al. 


isk, a nd passing stars was hrst done by 


Duncan et ah 


(120041) with more realistic initial conditions. 


Pones et al. 


(1987) 


(12004 


showed the time evolution of the mean eccentricity and inclination of planetesimals and 
their dependencies on the initial semimajor axes. In their calculation, the mean values 
expected for the isotropic distribution of planetesimals are attained in 4.5 Gyr, for the 
semimajor axis larger than 10^2x10^ AU. However, the time evolution of the distributions 
of the eccentricities and inclinations, and the role of stellar encounters in their evolution are 
not clarified. 


In the present paper we investigate the effect of stellar encounters on the evolution of a 
planetesimal disk, from the point of view of Oort Gloud formation. We obtain (i) the impact 
parameters that describe the evolution timescales of the orbital elements, (ii) the decay 
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timescales (e-folding time) of surviving comets in the solar system, and its radial variations, 


and (iii) the evolution of the semimajor axis distribution in the Oor 
classical impulse approximation (e.g.. 


Rickman 


1976 


Rickman et al. 


cloud. We use the 


20051 ) to calculate the 


velocity change of planetesimals induced by stellar encounters. The outline of this paper is 
as follows: We hrst describe the basic dynamics of planetesimals due to a stellar encounter 
in Section 2. In Section 3, we describe the numerical model of our simulations. The results 
are presented in Section 4. Section 5 is devoted to a summary and discussion. 
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2. BASIC DYNAMICS 


We use the cl assical impulse approximation to calculate the velocity change given by a 


stellar encounter flRickman 


19761) . In this approximation, a planetesimal is held hxed with 


respect to the Sun, a star passes with constant velocity along the straight line, and the 
time intervals before and after the closest approach are assumed to be inhnity. Then, the 
velocity change of the planetesimal is given by 


Au = 


2Gm* 


62 

p 


b 

IP 


( 1 ) 


where G is the gravitational constant, m* and u* are the stellar mass and velocity, and bp 
and b are the position vectors to the star from the planetesimal and the Sun, respectively. 


In this section, we derive the impact parameters that can be the indexes for the orbital 
evolution. We deal only with planetesimals that have the orbital velocity much lower than 
the typical passing star, i.e., the planetesimals with a > 10^ AU. 


2.1. Ejection 


Assuming bp b, we use the following simplihed equation instead of Equation ([T]), 

2Gm* 1 


Av ~ 


V* bp 


( 2 ) 


Substituting Av = Vesc = ■\/2GMQ/r into Equation ([2]), we obtain the impact 
parameter of the star-planetesimal encounter that gives the velocity change large enough to 
escape from the Solar system. 


m* / 2G 

“ “ iirvW’ 


(3) 


where r is the heliocentric distance of the planetesimal and we used v Av. This condition 
is basically satished when a planetesimal is near its aphelia. When bp 3> b, bp in Equation 
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(l2|) is substituted by b. Then we obtain the impact parameter of the star-Sun encounter that 
gives the velocity change large enough to escape from the Solar system to the planetesimal 
with r, which is identical to Equation (E]). This is the same to the other impact parameters 
derived in Sections 12.21 and 12.31 


Using besc, we derive the parameter dependence of the survival rate of planetesimals in 


the Solar system Pbound analytically by improving Weissman’s method fjWeissman 


1980|). 


Weissmaru f)l980l) estimated the fraction of the Oort cloud comets ejected by a single stellar 


encounter by assuming that a star drills a narrow tunnel through the Oort cloud, ejecting 
all the planetesimals within a radius of feesc- Assuming a constant number density of 
planetesimals in the Oort cloud and the length of the tunnel equal to the mean chord length, 


Weissmaru (jl980[) estimated that nine percent of the cloud population is ejected by a single 
encounter with a 1 Mq star with relative velocity 20 kms“^. Using the same assumption, 
we express the ratio of the number of planetesimals ejected by a stellar encounter to that of 
surviving planetesimals as 




7rP2 


oc a 


-1 


n* 


( 4 ) 


where a is the semimajor axis and we used r oc Pc oc a. The factor 2 of expresses 
the contributions from stellar encounters with planetesimals and the Sun. Additionally 
assuming that Rc is constant with time, the number of stellar encounters within 6 < Pc in 
time t is proportional to 


B oc /encPcf OC 


( 5 ) 


where /enc is fhe encounter frequency per 1 Myr within 1 pc from the Sun. Using A and P, 
Pbound is expressed as an exponential decay. 


fljomid = (1 - A)® ~ exp 



( 6 ) 
















where we used A <^1, and te is the e-folding time of the Oort cloud, 


te oc a 



( 7 ) 


2.2. Eccentricity change 


Next, we dehne the star-planetesimal impact parameter be, for which a planetesimal 
gains the eccentricity change Ae. We consider a planetesimal in the Cartesian coordinates 
{x, y, z) centered on the Sun and the x-axis is chosen as parallel to the position vector 
of the planetesimal from the Sun. The x-y plane agrees with its orbital plane (hereafter 
refered to as the reference plane). Then the position is written as r = (r, 0, 0). 

v^r — (v ■ r)v r 
^ ~ GMq r’ 

After a stellar encounter that gives the planetesimal the velocity change of An = 

{Avx, Ai.Vy, Avz), the planetesimal has the new eccentricity vector e'. The velocity change 
due to a stellar encounter is given by 



' Avx ' 


cos a ' 

An = 

AVy 

= Av 

sin a cos {3 


K An^ j 


y sin a sin {3 / 


(9) 


where cos a and (3 are randomly chosen from a uniform distribution between -1 and 1 and 0 
and 360°, respectively, since stellar encounters occur with random direction with respect to 
the Sun and the planetesimal. The change of the eccentricity vector is expressed as 


Ae 


e — e = 


GMq 


( {Avf - {Avx? \ 

-AvxAvy 

\ -AvxAvz / 


( 10 ) 
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where we used the approximation v Av. Then the absolute value of the eccentricity 
vector change is given as 


. rAv r/A n2 /a x2i1/2 r(Av)‘^ . 

From Equation flTT]) . we obtain the velocity change required for the change of Ae 


AuAe = 


Combining Equations (EJ and flT^ . we have 


GMq Ae 
r sin a 


( 11 ) 


( 12 ) 


h, = 2 


m* 


G r 


u* V Mq Ae 


sina ~ 


2 sin a 


(13) 


The ratio between the cross-sections derived from Equations ([2D and ffT2D is ~ 2 for 
sin a ~ Ae. This ratio implies that when one planetesimal escapes from the solar system 
due to a stellar encounter, another planetesimal reduces its eccentricity and evolves into 
a nearly circular orbit. Additionally, Equation (|T3ll tells us that Ae is Levy flight: the 
probability for a planetesimal having Ae > e by a single stellar enc ounter P(Ae> e), which 


is expected to be proportional to b^, follows P(Ae > e) oc e " (e.g.. 




Mandelbrot 


198211 . The 


other examples of the Levy flight are found in the semimajor axis evolution of planetesimals 


by planetary scattering (see Appendix. A) and in the evoluti on of the eccen 


inclination of the binary interacting with unbound perturbers fICollins Sz Sari 


ricity and 


20081 ). 


Next, we define the change of the argument of perihelion Au as 


cos Au) = 


e ■ e 


e e' 


(14) 


Under the assumptions of u ■ r = 0 and v Av, we have e ~ (—1, 0, 0) and then |e| ~ 1. 
Substituting v' = v + Av into Equation ([HD, we obtain 


/ 


e ~ 


iAv)^r 

GMq 

(Av)^r 


sin"^ a — 1 


\ 


r . n 

cos a sm a cos p 


■ GMq 


COS a sm a sm 


(15) 


ta/3 / 
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When Aco = 90°, 


, (Avfr . 2 „ 

e ■ e = 1 — „, , — sin a = 0 


(16) 


GMq 

Snbstitnting Eqnation (|2]) into Eqnation (1T6|) . we obtain the impact parameter 6^ that gives 
the planetesimal Au = 90°, 


= 2 - 


G 


' sin a = y/2 sin abe 


n* \l Mq 

Eqnation (fT7)) tells ns that h^j is on the same order of be for Ae = 1 bnt slightly smaller 
than it since sin a < 1 and the mean valne of sin a is vr/d < 1. 


(17) 


2.3. Inclination change 


h' = r X {v + Av) = 


cos Afl sin Ai 
cos Ai 


(18) 


We dehne Ai as an inclination of the planetesimal against the reference plane after 
a stellar enconnter. Under the same assnmption nsed in the derivation of be, the angnlar 
momentnm of the planetesimal after a stellar enconnter is written as 

^ 0 \ f sinAflsinAi ^ 

—rAvz = 

r{v + Avy) y y cosAi 

where Afl is the longitnde of the ascending node from the a:-axis. Under the impulse 
approximation, h'^, the a:-component of h' , is always 0 because the a:-axis lies on the orbital 
plane of the planetesimal. Using h'^ = h'sinAUsinAi = 0 in Equation fflSj) . we have 
Aff = 0, which means is not changed by a velocity kick given at the aphelion. Thus, Ai 
can be written as 


f ^' v \ ( ^'^2 

Ai = arctan -^ = arctan --— 

\h'J \v^AVy 

Assuming Av n in Equation flT^ . and using Equation ([9]), we have 

B (/3<180°) 

Ai= \ 

/3 - 180° {/3 > 180°) 


( 19 ) 
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where is uniformly distributed between 0 and 360°. 


Here we use Av = u as a condition to derive the impact parameter bi, which gives the 
planetesimal Ai = 90°. The orbital velocity of the planetesimal at its aphelion vq is written 
as 




' 2 GMo(l-e) 


Q 


( 20 ) 


where Q is the aphelion distance. Using Av = vq and r = Q in Equation ([2]), we obtain 


, rn* 

bi = 2 — ^ 
u* 


G 


( 21 ) 


Ms(l-e) y(l-e)"“” 

which means that for 1 — e ~ 10“^, when one planetesimal escapes from the solar system 
due to a stellar encounter, other ~ 2000 planetesimals gain ~ 90° inclinations. Comparing 
bi to fee or we hnd that the timescale the f-distribution requires to relax is about three 
orders of magnitude shorter than that for the e-distribution. 


3. METHOD OF CALCULATION 

We describe the initial distribution of planetesimals and stellar parameters, and the 
orbital evolution of planetesimals by stellar encounters. The parameter sets of all models 
are summarized in Tables [H 121 and El 


3.1. Planetesimal Disks 

We set up initial planetesimal disks formed by the transportation of planetesimals due 
to planetary scattering. In such planetesimal disks, planetesimals have their perihelion 
distances in the planetary region with small inclinations. We performed the test calculations 
and found that the disk evolution depends mainly on the semimajor axis, and weakly on the 
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perihelion distance if it is inside the planetary region and a > 10^ AU (i.e., e ~ 1 ). Thus, for 
simplicity, we consider planetesimals with perihelion distances go = 10 AU and inclinations 
io = 0 as the standard ( 10 ) model (hereafter we use the subscript 0 for the initial value). 
The initial semimajor axes of the planetesimals are oq = 5 x 10^, 10^, 1.5x10^, ..., and 
5 X 10^ AU. Their initial angle variables, the argument of perihelion wq, the longitude 
of ascending node flo) cind the mean anomaly Mq are uniformly distributed in the range 
0-360°. Additionally, we prepare spherical Oort clouds that have uniform distributions 
of — 1 < cosi < 1 and 10 AU< g < oq as initial conditions to evaluate the effect of the 
cloud shape on the decay time scale (12 in Table [T]). Each disk or cloud consists of 10^ 
planetesimals with identical oq. 


To investigate the radial variation of the structure and evolution of the whole disk, we 
consider planetesimal disks that have broad oo-distributions in 5 X 10^ AU< ao < 5 X 10“^ 
AU as summarized in Table The an-distributi on for the standard model (WO) follows 
dn/dao oc Oq, where 7 = —2 flHiguchi et al.ll2006f) . which is a flat distribution in E = 1/a 
(see Appendix). We test four additional disks for 7 = 0, 1, -3, and -1 for comparison (Wl, 
W2, W4, and W5). Each disk consists of 5 x 10® planetesimals. 


3.2. Passing Stars 


The stellar encounters are assumed to occur with random directions and follow the 
distribution dug/db oc 6 , where rig is a number of stellar encounters. The time interval 


to th e next stellar encounter is given according to a Poisson distribution flHeisler et ah 


19871 ). The star sets are described by hve parameters; the stellar mass m* and velocity n*, 
the encounter frequency /enc (in number per 1 Myr within 1 pc from the Sun), and the 
minimum and maximum impact parameters 6 min and &max- 
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We consider two types of star sets. First, we generate the star sets that consist of stars 
with identical m* chosen between 0.25Mq and 2Mq, and n* = 20 kms“^ (nnder the impnlse 
approximation, the important valne is and we can £x one of m* or n*). We choose 

m* = O.SMq as the standard model. The second star sets consist of stars with the realistic 
distributions of m* and n* with fenc for each type of star, based on the observations of the 
Solar neighborhood (13 in Table [Hand W4 in Table [2]). 


The stellar mass used in the standard model m* = O.SMq is estimated from the 
observations. Let us assume that the energy change due to a kick by fc-type star can be 
written as 


AEk oc (An)^ oc 




9{b, bp) 


2 


( 22 ) 


where we assume g{b,bp), the term given as a function of the impact parameters, is 
independent of k. The stellar velocity is given as 


'^*,k = '^0 + 3^^ “ 2\/3nQ(jcos6', 


(23) 


where Vq is the solar apex velocity, a is the ID velocity dispersion, and 0 is a direction 
angle between n* and a that has an uniform distribution between 0 and 360°. We assume 
that the typical parameters fenc, and h* satisfy the following equation, which is derived 
by averaging over 6 and summing up the energy changes from all types of stars, 

= fenc f ) • (24) 

k=0 ^ 

Adopting Vq, a, fenc and m* in Table [3l fenc = 10 and h*=20 kms“^, we have m* ~ O.SMq. 


To generate the distributions of m* and n*, we follow the procedure described in 


Rickman et ah 

(200§ 

) anc 

8 of 

Garcia-Sanchez et ah 


use t 


tie same va 


ues_of Uq and a , which are taken from Table 


(1200ll) and lAllenI fllQSSi) . This realistic star set consists of 13 


types of stars. The star data in this paper is summarized in Table [3l The sum of the 



























14 


encounter frequency of all 


19871: 


Rickman et al. 


2994 


types of stars is estimated to be /total — 10.5 (e.g., 


Heisler et al. 


29981) . Then one can estimate that one star approaches the Sun 


within 9.995 pc (~ 19^ AU) in 5 Gyr. Although we calculate the evolution for 19 Gyr, we 
set femin = 9.995 pc for all star sets and 6max = 1 pc for the standard model (19). It has been 
said that the eff ect of the stars for 6 > 1 pc is smaller than that of stars with 6 <1 pc (e.g., 


Weissman 


19891) . We also perform several calculations to evaluate the effect of distant stars 


with 6inax =9.25, 9.5 and 2 pc (models 14, 15, and 16) and hnd that the differences among 
&max = 9.25, 9.5,1, and 2 pc are small if 2ao < &max. Therefore, 6max = 1 pc is appropriate 
for the Oort cloud with = 19® AU. We can also show it analytically (see Appendix). 


We perform 19 runs for each model with different orientations and time of stellar 
passages that are randomly chosen. The error bars shown in the following hgures indicate 
the 1-a error from the 19 runs. 


3.3. Orbital Evolution 

We describe the procedures to calculate the orbital evolution due to stellar encounters 
and the Galactic tide. We develop a hybrid code that takes into account both stellar 
perturbations and the Galactic tide. 


3.3.1. Stellar encounters 

To calculate the orbital change of planetesimals due to the perturbation from a passing 
star, we use the classical impulse approximation that gives the velocity change described 
in Equation ([1]). For the impulse approximation to giye a reasonable result, the Kepler 
period of the planetesimal must be longer than the typical encounter time of a star with 
the Sun (~ 19® yr). We choose the minimum value of ao=5xl9® AU from this restriction. 
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Rickman et al 


(120051) showed that the classical impulse approximation gives a reasonably 


good approximation for the Oort cloud simulation as long as we avoid too close and/or too 
slow encounters. 


We calculate the orbital evolution of planetesimals by adding the velocity change for 
10 Gyr. If a planetesimal has its heliocentric distance r > 1 pc, it is counted as an escaper 
from the Solar system and discarded from the calculation. This criterion is not equivalent 
to e > 1, but we found that with some extra calculations that it is almost the same. The 
calculation is stopped when the survival rate Pbound = n/no < 10“^ or at 10 Gyr. 


3.3.2. Galactic tide 


We add the effect of the vertical component of the Galactic tidal force to the standard 
model (II in Table D and W3 in Table 2]) . This effect i s analy tically computed using the 


formulae derived in 


Higuchi et al. 


(120071) 


Higuchi et al. 


(120071) investigated the effect of the 


Galactic tide on the orbits of Oort cloud comets and obtained the evolution of the e, i, co, 
and G in the Galactic coordinates (note a does not change since the energy is conserved). 
They neglected the radial component of the Galactic tide so the energy and 2 ;-component of 
the angular momentum of comets are conservative quantities. The period of the oscillations 
of e, i, and u and the mean period of the circulation of G are as follows: 




vr 




n„ 


tn* — 


TT 


(25) 

(26) 


A3ll[K{k), ’ 

where is a mean motion of u* that is an angle variable directly related to uj, H 3 , a, 
and k are constants given by initial orbital elements of each planetesimal, H is an ecliptic 


integral of the third kind, and K{k) is a complete ecli ptic integral o 


derivation of the timescales and constants are given in 


Higuchi et al. 


the hrst kind. The 


( 2007 ). Note that in 
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Te,i,uj is written as P^j *, the factor of Equation (A24) is not correct Q, 
and t^*, which is written as Pq*, is given in a different form using a Fourier series. We 
assume that the timescale of the evolution due to the Galactic tide is oc or Pq. Both 

Pe,i,uj and Pn are proportional to n~}, which is 

Z/2 TT 

where uq is a vertical frequency of the solar motion in the Galactic disk that is constant, n' 
is the mean motion and Oq and 02 are the constants given by the initial orbital elements of 
the planetesimal. Then the timescale is proportional to The evolution timescales of 

these orbital elements except G are coupled. That means the Galactic tide itself does not 
randomize the planetesimal distribution. Though the evolution of G is independent of those 
of the other orbital elements, their periods are nearly commensurable. Therefore, at the 
beginning of the evolution, i and G are nearly coupled and make two strong peaks in the 
inclination distribution around 27° and 153° with respect to the ecliptic plane, where the 
inclination of the ecliptic plane to the Galactic disk is assumed to be 63°. 


Higuchi et al. 


(2007h 


In the simulations, planetesimals are affected by the Galactic tide during the time 
intervals between one stellar encounter and the next. After each stellar encounter, the new 
conservative quantities of planetesi mals are recalculated. We adopt the total density in the 


solar neighborhood p = O.IMqPC flHolmberg fc Flvnn 


2000 ). 


^ The factor A 2 in Equation (A24) in 

2 


Higuchi et al. 


Aq — 


zAo — 


3 Xq -3 


( 20071) should be A 3 , 
1 


Ai\/a2 ~ 0^0 2 Po ~ P \J0t2 — olq 


(27) 
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4. RESULTS 

4.1. Planetesimal Disks with Identical oq 

We present the results of simulations for initial planetesimal disks that consist of 
planetesimals with identical Oq. First, we show the result for the standard model with 
oq = 2 X 10“^AU. Next, we compare the results for different Oq and those with the Galactic 
tidal force. 


4 . 1 . 1 . Evolution of Distributions 

Figure m shows snapshots of the planetesimal disk plotted on the x-z plane of the 
Cartesian coordinates (x, y, z) for t = 1 Myr, 5 Myr, 10 Myr, ..., and 10 Gyr. The 
planetesimal disk is initially on the z = 0 plane. The planetesimals initially have oq = 2x10^ 
AU. By 10 Myr, the disk is almost flat. Some asymmetric structures seen at 5 Myr and 10 
Myr are the result of recent stellar encounters. At 100 Myr, the disk slightly expands from 
the initial size and has thickness of ~ 10^ AU. The planetesimal distribution expands to 
around 10^ AU and becomes almost spherically symmetric at 5 Gyr and shrinks by 10 Gyr. 

For quantitative analysis of the structure evolution, we introduce an indicator Rp that 
is the heliocentric distance that contains p percent of surviving planetesimals inside it. The 
time evolution of Rgo scaled by oq (i.e., i?ioo — 2ao at f = 0) is plotted in Figure El The 
structure becomes the largest in Rgo at ~ 2 Gyr for oq = 2 x 10^ AU. After reaching the 
maximum, it shrinks. The time evolution of Rgo for the other parameters is also plotted 
in Figure El The evolution of Rgo for m* = IMq is quantitatively the same as that for 
m* = O.SMq, but it is much more rapid, reaching the maximum value at ~1 Gyr. The 
evolution for m* = 0.25 Mq is slow and does not show the decrease in Rgo within 10 Gyr. 
The dependence of the maximum Rgo on oq is not well scaled by qq. The larger oo has the 
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shorter evolution timescale and smaller maximum value of Rgo/ao. This is simply because 
the maximum heliocentric distance for bound planetesimals is r^sc = 1 pc in our model, 
which is independent of Qq. 

Figures [3}l6] show the time evolution of distributions of semimajor axes, eccentricities, 
inclinations, and arguments of perihelion of the planetesimals plotted in Figure [H 
respectively. In 5 Gyr, the semimajor axis distribution evolves approximately into 
dn/da oc a“^. The diffusion is not symmetric about the initial value. At the beginning of 
the evolution, there are more planetesimals for a > ag. The mean value of a shifts and it 
reaches ~ 4 x 10^ AU. After the outward shift, the a- distribution finally shifts inward and 
the mean value of a decreases to ~ 3 x 10^ AU at 10 Gyr, since planetesimals with larger a 
are easily removed from the disk. 

The short-dash curve plotted in the e- and i- distributions is the distribution expected 
for an isotropic distribution. The eccentricity distribution in Figure H] gradually relaxes 
and attains the isotropic distribution, which is proportional to e. On the other hand, the 
inclination distribution shown in Figure 0 relaxes much more rapidly compared to those of 
the other orbital elements. Only in 5 Myr, the distribution already has its range from 0 to 
180°. This rapid spread of ^-distribution compared to the other orbital elements and the 
flat distribution are explained in Equation fl20l) in section 12.31 They are at the stage where 
only the i- distribution has relaxed its initial distribution, but the other orbital element 
distributions have not yet. The i- distribution ranging from 0 to 180° does not mean the 
Oort cloud is spherical. To be a spherical Oort cloud, the isotropic distribution of oj is 
also required. Figure [6] shows the time evolution of the cn—distribution. In 5 Myr, the 
(u-distribution is concentrated at uj =0 (or oj = 180° if i > 90°) i.e., the eccentricity vectors 
are still close to the initial plane. 

Figure [7] compares the evolution of ^-distribution due to the passing stars only (10, 
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left), the Galactic tide only (GT, center), and both of them (II, right) for oq = 2 x 10^ AU 
at 50 Myr, 500 Myr, and 5 Gyr. Two peaks in the ^-distribution seen in GT (center) are the 
strong feature made by the Galactic tide. They are also seen in 500 Myr of the II model. 
The distribution for II is just intermediate between those with and without the Galactic 
tide. At 5 Gyr, the difference between the distributions with and without the Galactic 
tide is quite small. The random-walk nature of passing stars can randomize the angular 
momenta of planetesimals quickly, with much shorter timescale than that for ejection, as 
shown in section [231 


4.1.2. Decay of Comet Clouds 


As described in the previous section, the number of planetesimals decreases while the e- 
and ^-distributions approach the isotropic distribution. We £t the decay curve empirically 
using the standard expo nential decay curve and t he stretched exponential decay defined by 


the Kohlrausch formula (IDobrovolskis et ah 


2nn7i . 


The empirical £t for Pbound obtained by the least-square £t optimized for uq = 2 x 10"^ 
AU is 


P 


"bound exp 


t 


(29) 


~ 5.6 


Oq 


2 X 10^ AU 


- 1.4 


m* 20 kms 
0.5Mq u* 


- 1 \ - 1-7 


/ei 


10 Myr 


-1 


Gyr. 


(30) 


We compar e Equ ation (l30|) to a previous study on the lifetime of wide binaries. 


Weinberg et al 


fjl987l ) analytically and numerically calculated the lifetime of a binary as a 


function of the o rigina l binary semimajor axis (oq) and other parameters of the perturbers. 


Weinberg et al. 


fll987l) showed the characteristic lifetim e of a bina r y wit h relatively small 


initial semimajor axis (< 0.1 pc) is proportional to a. 


-1 

0 ■ 


Tremaine 


(119931 ) referred to this to 
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describe the lifetime of the Oort cloud. 


Weinberg et ah 


fjl987() also show that the lifetime 


is proportional to where n* is the number density of perturbers (oc fenc/v*), 

M* is the mass of a perturber, and l^ei is the relative velocity between the binary and the 
perturber. The dependence on oq in Equatio n fl30ll is close to -1.34 , which is for a wide 


binary with the separation > 0.1 pc shown in 


Weinberg et ah 


(119871). 


The curves in Figure [8] are Equation fl29|) using Equation fl30l) for each parameter set. 
They agree well with Pbound- The differences between Equations ([7]) and fl30|l are basically 
due to the two assumptions for Equation ([71) : (1) Pc is constant over the evolution, and (2) 
Av S> u in derivation of 6esc- The reason the dependence on oq in Equation fl5U]l is stronger 
than that in Equation (I7|) may be due to the evolution of a. The assumption of the constant 
Pc(oc Oq) is obviously broken, as seen in the evolution of Rqq in Figure 121 As a evolves, the 
evolution of the Oort cloud is accelerated and the dependence on oq becomes stronger. 


The stretched exponential decay is given by the Kohlrausch formula. 


Pbound exp 


/3' 


(31) 


where tn is a cons 


flDobrovolskis et af 


ant a nd (3 is the index known as the stretching parameter 


20071) . When [3 = 1, Equation fl3Tl) is the standard exponential 


decay and R = te- The e-folding time is given as 

Pbound 


te = 




(32) 


|dPbound/dt| 

which indicates that decreases and increases with time for [3 > 1 and [3 <1, respectively. 
We calculate R and (3 using the least-square £t and obtain 

m* 20kms“^\ / /e. 


fJl* = 5.3 


Oq 

2 X 10^ AU 


- 1.4 


0.5Mq u* 


10 Myr 


-1 


- 0.97 


Gyr, (33) 


and 




2.0PA^ 


bound’ 


(34) 
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which shows that Equation fl3T]) is now an implicit function of Pbound- As seen in Equation 
flM|l . (3 is not a constant over the evolution. Initially (3 is much larger than 1 and decreases 
with Pbound and when Pbound — 0.18, /? ~ 1. The fits expressed by Equations fl5^ and fIMD 
agree much better than those with the standard exponential decay, Equation fl30|) . 

All the analytical arguments in Section [2] are done under the assumption that the 
Oort cloud is spherical, despite the initial flat distributions of planetesimals. To see how it 
affects the results, we also perform the same calculations for the initially spherical cloud 
(12 in Table [I]). The e-folding time of the 12 model is plotted against oq in Figure |9l We 
found that the e-folding times of the 12 and 10 models are almost indistinguishable and the 
effect of the initial flat structure is negligible in Pbound under the random and many stellar 
encounters. 

Figure [To] shows Pbound for the realistic star sets (13) against time with Pbound f^^ 

/enc = 10, u*=20 kms“^, and m* = O.SM©. The error bars are large due to stochastic close 
encounters of massive stars, mostly BO stars. However, the estimation of the evolution 
due to the realistic stars using P^o^jd for / = /encW* = h*, and m* = m* agrees with 
Pbound within l-cr. Figure [TTl shows the time te against Oq with 1 -a error bars and for 
m* = O.SMq. They also agree within l-cr except for uq = lO'^ AU. This agreement between 
10 and 13 shows that the averaging of energy kicks from many types of stars described in 
Section 13.21 is a good approximation. 

The Galactic tide makes little difference in tg as seen in Figure [9l It is easy to 
understand, because the Galactic tide in our model does not change the energy of 
planetesimals or the cross-section of the Oort cloud (~ Qq). The timescale of the evolution 
due to the Galactic tide is proportional to (Eq. (I28|) b This power-law index is very 

close to that of the evolution due to stellar encounters in Equation fl30l) . which is -1.4. 
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4.2. Planetesimal Disks with ao-Distribution 

We show the evolution and lifetime of disks that consist of planetesimals with 
ao-distribution and several additional effects. 


4 . 2 . 1 . Decay of Comet Clouds 

Figure [12] shows Pbound for the disks with 7 = —2 and 0 (WO and Wl). The e-folding 
time obtained by htting are 15.9 Gyr and 5.2 Gyr for 7 = —2 and 0, respectively. We 
approximate the decay of the planetesimal number for the disk with oo-distribution 
using inbound planetesimal disk of identical Qq. The decay curve for the disk with 

ao-distribution is obtained by averaging P^ound or ([3TD) over oq but it cannot 

be obtained analytically. Therefore, we estimate the decay curve by averaging Equation 
fl30|l instead of Equation fl2^ . We assume that planetesimals have ao-distribution following 
dn/dao cc Uq in the range of 5 x 10^ — 5 x 10"^ AU. Integrating Equation fl30|) over ao using 
the probability density distribution of ao that is proportional to dn/dao, we obtain the 
averaged e-folding time (te) 17.9 Gyr, 6.7 Gyr, and 3.9 Gyr for 7 = —2, 0, and 1 (WO, 
Wl, and W2. See Table [2D, respectively. They are slightly shorter than the e-folding times 
obtained by htting. The agreements of Pbound and Abound shown in Figure [12] are reasonably 
good as well as those for the identical ao models (Fig. [HD- 


4 . 2 . 2 . Evolution of the disk structure 

Figure [13] shows the time evolution of r-distributions for the disks WO and Wl, 
respectively. The disk WO keeps the initial distribution during the evolution. On the other 
hand, the r-distribution for the disk Wl greatly changes and approaches that with 7 ~ —2 
where we assume that the disks have the r-distributions that follow dn/dr oc r"^. Figure [13] 


23 


shows the r-distributions in 10 Gyr for a 


are close to dn/dr oc r 


-2 


Rickman et ah 


1 the d isk models. All the distribution in 10 Gyr 


( 120081 ) also found in the calculation that the disk 


initially with 7 = —1.5 evolves into a distribution with 7 ~ —2. The final r-distributions 
with 7 = —2 can be explained as follows: The energy kicks that are small but repeatedly 
given by stars, relax the initial energy distribution and result in a wide and smooth energy 
distribution whose width is larger than the energy range of the Oort cloud. Therefore, when 

he Oort cloud, the energy dis tribution is roughly 


applied to the narrow energy range of 
flat, from which follows dn/da oc flMalvshkin &: Tremaine 


19991). Using a oc r, we 


expect r-distribution follows dn/dr oc r 


-2 


Next, we divide the disk evenly by r and see how the disk evolves into a spherical 
structure from a flat disk. To evaluate this structural evolution, we introduce a new index 
which is a normalized ratio of the sums of the square of the vertical axis and the radial 
axis defined by 


S + 1/2) ’ 


(35) 


where x,y,z are the coordinates that satisfy < x^ + < {r + Ar)^. When the 

structure between r and r -|- Ar is perfectly flat, 0 ^= 0 . The r-dependence of for WO 
is plotted in Figure [H] at various times, showing the evolution of the flat disk into the 
spherical cloud. Since the flare-up of the disk occurs contiguously with r and time, it is 
difficult to tell from this hgure where the boundary of the inner (relatively flat) and outer 
(spherical and isotropic) Oort clouds is. The high at r < 5 x 10^ AU is due to new 
comets injected from the outer Oort cloud. 


Figure HH shows the r-dependence of for the model with the Galactic tide (W3). 
We can see that the Galactic tide (1) accelerates the flare-up of the disk, especially at the 
beginning, and (2) produces a wave-like structure propagating inward. The wave is due to 
the nature of the Galactic tide that causes a periodic oscillation of the inclinations against 
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the Galactic plane {i = 63°). Therefore all planetesimals follow similar orbital evolution 
with different timescales that depend on the orbital elements (especially on a ). The 
effect of passing stars is too weak to let planetesimals forget their conservative quantities 
completely within the timescale of the oscillation due to the Galactic tide. However, their 
inclinations do not have an isotropic distribution even for a = 0.5 since the Galactic tide 
does not randomize the orbital elements of planetesimals. Due to the wave, it is difficult to 
distinguish the inner and outer Oort clouds using in W3. 


4-2.3. Radial transport of planetesimals 

We calculate the decay of planetesimals in each r range Pbound(?") and derive their fits. 
However, due to the radial transport of planetesimals, not all Pbound(^") decay monotonically. 
We deal with only r-bins whose Tbound(^") shows a monotonic decay. For WO, Tbound(^) only 
r < 8 X 10^AU shows a monotonic decay, while Pbound for r > 8 x 10^AU region shows an 
increase in the number of planetesimals due to the outward transport of planetesimals. 
Assuming that the decay curves follow the standard exponential and Kohlrausch formula, 
we calculate te, to, and ft for r < lO^AU by fitting. Also for the r-bins that do not have a 
monotonic decay, we obtain t^, to, and (4 using the decreasing part of Pbound- 

First, we compare the results for disks with different 7 (models WO and Wl). Figure ?? 
shows averaged fg and to obtained by the fitting of numerical results with error bars against 
r. For 7 = —2, tg and to are almost flat at tg ~ to ~ 10 Gyr for any r. In contrast, tg and 
to for 7 = 0 decrease inversely proportionally to r. This indicates that the efficiency of the 
planetesimal supply from the inner Oort cloud to the outer Oort cloud strongly depends on 
the oo-distribution. The dependence of tg oc r~^ corresponds to Equation ([7]) derived under 
an assumption that a-distribution is constant (i.e., no planetesimal supply). Figure ?? also 
shows ft obtained by fitting. We find that WO and Wl have /3 > 1 for the inner Oort cloud 
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and /5 < 1 for the outer Oort cloud fsee I4.1.2|) . This means that no matter how massive the 
outer Oort cloud is, the net rate of transportation of comets due to the perturbation from 
stars is outward. 


Second, we compare the evolution for the disk WO to those under the Galactic tide 
and non-identical stars (W3 and W4) in Figure ??. In general, the effects of the Galactic 
tide and the non-identical stars are small although the error bars in model W4 are much 
larger than those for the other disks. But in detail, we can see that the Galactic tide 
increases te and to and reduces (3. This slight stabilization may be due to the decrease of 
mean heliocentric distances of planetesimals, a(l -|- e^/2), as a result of the decrease of e of 
planetesimals that are in the Kozai-Lidov cycle. In the Ko zai-Lidov cycle, the eccentricity 


of most of the planetesimals initially with e ~ 1 decreases flHiguchi et ah 


200711 . 
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5. SUMMARY AND DISCUSSION 

We have investigated the effect of stellar encounters on Oort cloud formation from 
a flat disk to a spherical structure. First we analytically derived the impact parameters 
for e, and u changes and compared them to the impact parameter for ejection using 
the classical impulse approximation. Using these impact parameters, we showed that the 
relaxation timescale for i is much shorter than those for e, u and ejection and e does not 
evolve according to a normal random walk but a Levy flight. 

Next, we performed numerical calculations of the evolution of planetesimal disks 
by stellar encounters using the classical impulse approximation. We numerically and 
analytically showed that the encounter with distant stars is not effective, either in ejection 
of planetesimals or in randomization of the orbital elements. We empirically derived 
the htting formulae Rbound using the standard exponential formula and the Kohlrausch 
formula. We obtained the timescales te and to uud the stretching parameter /? for the disk 
that consists of planetesimals with identical initial semimajor axis Oq and the local te, to, 
and (3 as functions of the heliocentric distance r. Using te, to, and (3, we discussed the 
timescale and efficiency of transporting of comets from the inner/outer Oort cloud. We 
found that the timescale of the decay of the Oort cloud is almost the same at any r, as 
seen in the r-distribution that follows dn/dr oc during evolution. This is due to the 
sufficient supply of planetesimals from the inner Oort cloud to the outer Oort cloud. If 
the initial planetesimal disk has a relatively outer-massive distribution and the outward 
supply of planetesimals is small, the outer Oort cloud has a much shorter lifetime roughly 
proportional to r~^. We found that, in any disk model, the decay rate of the outer Oort 
cloud is smaller than the standard exponential decay rate (/9 < 1) because of the sufficient 
supply of planetesimals from the inner Oort cloud. We showed that the distribution of the 
semimajor axes of comets in the Oort cloud approximately follows dn/da oc for any 
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initial distribution of the semimajor axis. This corresponds to the flat distribution of orbital 
energies. Since the initial a-distribution of the sca ttered planete s imal disk generated by 


planet scattering follows dn/da oc a ^ as shown in 


Higuchi et ah 


teOOTI) . the a-distribution 


of the Oort cloud does not change during evolution. 


The local Galactic environment in the stage of Oort cloud formation such as stellar 



embedded in a giant molecular cloud just after it formed, the high density of the cluster 
accelerates the evolution of the Oort cloud. Investigation of the time evolution of the 
Galactic environment is left for future work. 


The distribution of the orbital elements, especially i of new comets produced by the 
Galactic tide is an interesting subject that is directly connected to observations. Our 
study showed that the timescales of the evolution due to stellar encounters is analytically 
proportional to (Eq. [7]) and numerically proportional to (Eq. [30] and independent 
of r (Fig. ??). As the timescale of the evolution due to the Galactic tide is proportional to 
(2-3/2 |28|) . the difference of the effects of stellar encounters and the Galactic tide might 

be seen especially in the anisotropic i- distribution of the new comets against their original 
semimajor axes. The production and the distribution of new comets will be the subject of 
our next work. 
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Initial a-Distribution of Planetesimals 

numerically showed that the probability distribution function for 
for the Oort cloud comet candidates formed by planetary scattering (i.e., widely distributed 
scattered planetesimal disk) is 

P{a > Oc) oc a“^, (Al) 


Higuchi et al. 


fl2nn6h 


where Oc is the minimum semimajor axis for the Oort cloud comet candidates. This is valid 
for ttc > lOopianet, wliere Opianet IS the semimajor axis of the planet that formed the Oort 
cloud comet candidates by planetesimal scattering. Thi s function means the scattering 


process of planetesimals by a planet is a Levy flight (e.g., 


Mandelbrot 


1982 ). Equation flAlj) 


shows that the distribution of the semimajor axes follows, 

dn 


oc a 


-2 


da 


which corresponds to the energy distribution 


dn 

— = const., 


(A2) 


(A3) 


where E = 1/a. This uniform distribution can be explained using a Gaussian energy kick 
distribution. If the distribution of the energy kick given by a planet follows a Gaussian 
distribution the energy distribution of planetesimals is 

Ht? 

— oc exp [-{E - Eof] , (A4) 


where Eq = 1 /ao and oq is the semimajor axis of the planetesimal before planetary 
scattering. The energy range for Oort cloud comets is very narrow and close to E = 0, 
which means the Oort cloud energy range lies on the low-energy tail of the Gaussian energy 
distribution. The gradient of the distribution is given by 

(^) ^ ^ ~ ~ ’ 
which goes to zero when Eq^ E and E —)■ 0. This means that 7 = —2 is expected. 


(AS) 
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B. Effect of Distant Stars 

We analytically derive the effect of passing stars with fep > besc, which does not eject 
planetesimals by a single enconnter. Assnming that the initial velocity of a planetesimal 
V is mnch smaller than An, the energy change per nnit mass dne to a stellar enconnter is 
given as AE ~ i(An)^. Then the energy change given by stars for bmin < b < bmux per nnit 
time is 

rfemax 

AEt = / 2TTbfAEdb, (Bl) 

^ ^min 

where / is a flnx of stars per nnit area per nnit time (i.e., / = /enc7i'~^Myr“^). Assnming 
bp <^b in Eqnation ([1]), we have AE oc b~'^. This means that AEt converges with b —)■ oo. 

To compare the effects of stars that do not penetrate the Oort clond and that penetrate 
it (bnt not close enongh to eject a comet by a single enconnter), we evalnate these effects 
as a fnnction of r. The averaged energy change per nnit time given by random stellar 
enconnters for fcgsc < & < oo is given by 

1 pi p7T poo 

(AEt) = — / / 27rbfAEdbd(3dcosa, (B2) 

^ Jo Jo Jbesc 

where a and /3 are the direction angles of stellar enconnters and cos a and /3 have nniform 
distribntions. We integrate Eqnation flB2j) neglecting the small terms and avoiding the 
region where a single stellar enconnter ejects comets (see Appendix. C). Then we have 


A£, = 0 = 1,2,3) 

(B3) 

0 i~ln? 7 , [b<r-besc\ 

(B4) 

©2 ~ A [r ^esc <b<r + &esc] 

(B5) 

©3 ~ 1 - 7r\/277"3, [r + besc<b] 

(B6) 


where rj = r/{m^/v^YMQ/2G = T^/feesc ~ O.SMq and n* = 20kms r/ S> 1. 

Assnming r as the radins of the Oort clond, we can say that Eqnations flB4p and fIBSp are for 
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stellar encounters penetrating the Oort cloud {b < r + besc), and Equation flB6j) is for that 
not penetrating (6 > r + 6esc)- We also numerically integrate Equation flB2p for femin = &esc 
and 6max=l pc and hnd that Equations flB4p - flB6p are good approximations within ~20% 
error for r > 10^ AU. Since the ratio (0i + 02)/03 > 10 for any r for m* = O.SMq and 
n* = 20 kms“^ (i.e., r/ ~ r), we conclude that the effect of stars not penetrating the Oort 
cloud is small. 


C. Integration of Equation (IB2p 


To perform the integration of Equation flB2H . we introduce coordinates similar to iBailv 


fjl983h : centered on the Sun with ; 2 -axis parallel to n* and choose the x-axis so that r is on 


the x-z plane. The last axis in the left-handed coordinate system is the j/-axis. Then h is 
on the x-y plane. Let the angle between r and j/-axis a and the angle between h and x-axis 
/5. When the direction angle of stellar encounters is isotropic, cos a and jS have uniform 
distributions between -1 and 1 and 0 and 360°, respectively. The stellar impact parameters 
against the Sun b and a planetesimal bp are 

^ bcos/3 \ 

6 sin /3 
0 

b cos (3 — r sm a 
6 sin (3 
0 


V 


(Cl) 


/ 


/ 


K = 


\ 


\ 


(C2) 


Equation ([T]) is rewritten as 


An = 


2Gm* 




/ bcos/3—rsma b cos fi \ 
6 sin (3 bsin/3 


~1W 


'W 


Vo 


(C3) 
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where 


sin^ a — 2hr sin a cos (3. 


(C4) 


The averaged energy change per unit times given by stellar encounters are given by 
Equation flB2p . We avoid the range bp < 6esc from the integral where a single energy kick is 
large enough to eject a planetesimal. Then b is divided into three ranges: (l)[6esc, r — 6esc], 
(2)[r — besc, r + 6esc], and (3)[r + 6esc, cxo]. We integrate Equation (1B2p by neglecting the 
{besc/bY terms and using 1 and obtained the following formulae for each b range, 


AE, = 


-©i (J = 1,2,3) 


(C5) 


01 ~ ^2\/2 — Inr^J r; 4 - r] 

~ Inr; 


+ In 


4^3/2 f 1 — '/ 2 ri 4 


(r^-1)1/2 

[b <r - be: 
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2 + - 
TT 


(C6) 


00 ~ 


2\/2 _5 ^ _3 \/2(2 + 7r) _1 , l + r7-^/2 

-p 4 + 2 '/ 2 ri 4 -I- rj 4 -)- 2 In 


TT 

\/2(2 + 7r; 

TT 


TT 




1 — r7“i/2 

[r - besc <b<r + 6esc] 


(C7) 


03 ^ 


1 — y^(l + ~ 1 arcsin [(1 + rj ^'^^) 


TT 




[r + besc < b] 


(C8) 
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Fig. 2.— Disk radius i?9o% that contains 90 percent of planetesimals inside it against t for 
model 10 of oq = 2 X 10"^ AU and m* = 0.5 Mq (solid), oq = 2 x 10^ AU and m* = 0.5 Mq 
(dashed), oq = 2 x 10^ AU and m* = 0.5 Mq (short-dashed), oq = 2 x 10^ AU and 
m* = 0.25 Mq (dotted), and Oq = 2 x 10^ AU and m* = 1 Mq (dot-dashed). 
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Fig. 3.— Normalized distributions of semimajor axes for model 10 of oq = 2 x 10^ AU at 
t = 1 Myr, 5 Myr, 10 Myr, 50 Myr, ... , and 10 Gyr. 
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Fig. 4.— Same as Figj3] but for eccentricity. The short-dashed lines show the isotropic 
distribution. 
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Fig. 5.— Same as Figj3] but for inclination. The short-dashed curves show the isotropic 
distribution. 
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Fig. 6.— Same as Figj3]but for argument of perihelion 
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Fig. 7.— Evolution of inclination distribution due to passing stars only (left, model 10), the 
Galactic tide only (middle), and both of them (right, model II) at 50 Myr, 500 Myr, and 5 
Gyr for oq = 2 x 10^ AU. 
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Fig. 8.— Surviving rate of planetesimals Pbound with 1-cr error bars against t for model 10 
of oq = 10^ AU and m* = 0.5 Mq (crosses), oq = 2 x 10"^ AU and m* = 0.5 M© (circles), 
Oq = 4 X 10"^ AU and m* = 0.5 Mq (open circles), Oq = 2 x lO'^ AU and m* = 0.25 Mq (open 
squares), and Oq = 2 x lO'^ AU and m* = 1 Mq (squares) with using the standard 

exponential decay (S, double short-dashed curve) and the stretched exponential decay (K, 
dotted curve). 
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Fig. 9.— e-folding time with 1-a error bars against Oq for models 10 (crosses), II (circles), 
12 (open circles), 14 (open sqnares), 15 (sqnares), and 16 (open triangles) with (donble 
short-dashed line). 
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Fig. 10.— Surviving rate of planetesimals Pbound due to encounters with non-identical mass 
stars (model 13) with 1-cr error bars for model 13 of Oq = 10"^ AU and m* = 0.5 Mq (crosses), 
oo = 2 X 10"^ AU and m* = 0.5 Mq (circles), and oq = 4 x 10"^ AU and m* = 0.5 Mq 
(open circles) against t with inbound '^* — 0.5Mq using the standard exponential decay 
(S, double short-dashed curve) and the stretched exponential decay (K, dotted curve). 
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Fig. 11.— e-folding time against Oq with 1-a error bars for the non-identical mass star 
model (model 13) and for m* = 0.5Mq (double short-dashed line). 
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Fig. 12.— Surviving rate of planetesimals Pbound against t for models WO (crosses), W1 
(circles), and W2 (open circles) with 1-cr error bars and Abound using (te) for WO (solid 
curve), W1 (dashed curve), and W2 (short-dashed curve). 
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Fig. 13.— Time evolution of distributions of heliocentric distances for models WO (top) and 
W1 (middle). The vertical axis shows the number of planetesimals in each r-bin. Distri¬ 
butions of heliocentric distances at f = 10 Gyr for all the models summarized in Table [2] 
(bottom). 
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Fig. 14.— Normalized ratio of the vertical and radial axes against r with 1-a error bars 
for models WO (top) and W4 (bottom) at t=100 Myr (crosses), 500 Myr (circles), 1 Gyr 
(open circles), and 5 Gyr (open squares). 
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Fig. 15.— Standard e-folding time against r (left), the time to for the Kohlrausch formula 
against r (middle), and the stretching parameter P for the Kohlrausch formula against r 
(right) for models WO, Wl, W3, and W4 from the top. The data without hlled circles are 
for the r-bins that do not have a monotonic decay. 










Table 1: Model Parameters for Identical Oq Disks 


Model 

ao [AU] 

go [AU] 

*0 [deg] 

m* [Mq] /enc [Myr 


&max [pc] 

GT 

10 

2 

X 

10^ 

10 

0 

0.5 

10 

1 

no 

(variations) 

5 X 10^5 

X 

10^ 

10 

0 

0.25-2 2- 

-20 

1 

no 

11 

2 

X 

10^ 

10 

0 

0.5 

10 

1 

yes 

12 

2 

X 

10^ 

10-2 X 10^ 

0-180 

0.5 

10 

1 

no 

13 

2 

X 

10^ 

10 

0 

Table [3] 


1 

no 

14 

2 

X 

10^ 

10 

0 

0.5 

10 

0.25 

no 

15 

2 

X 

10^ 

10 

0 

0.5 

10 

0.5 

no 

16 

2 

X 

10^ 

10 

0 

0.5 

10 

2 

no 


Note. — For all, no = 10"^, = 20 kms 6min = 0.005pc, and 6max = Ipc. The *o-distribution in 12 


follows dn/d*o oc cosjq (isotropic). GT means the additional Galactic tidal force. 
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Table 2: Model Parameters for Disks with ao-distribution 


Model 

7 

m* [Mq] /enc [Myr 

GT 

WO 

-2 

0.5 

10 

no 

W1 

0 

0.5 

10 

no 

W2 

1 

0.5 

10 

no 

W3 

-2 

0.5 

10 

yes 

W4 

-2 

Table [3] 


no 

W5 

-3 

0.5 

10 

no 

W6 

-1 

0.5 

10 

no 


Note. — For all, no = 5x 10®, qq = 5x 10® —5x lO'^AU, qo = lOAU, io = 0, = 20 kms 6niin = 0.005pc, 


and 5niax = Ipc. GT means the additional Galactic tidal force. 
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Table 3: Stellar parameters made from Table 1 in 


Rickman et al 


(1200811 ■ 


Stellar type 

k 

Vq [kms 1] 

a [kms 1] 

/enc [Myr"^] 

m* [Mq] 

BO 

1 

18.6 

8.5 

0.005 

9 

AO 

2 

17.1 

11.4 

0.03 

3.2 

A5 

3 

13.7 

13.7 

0.04 

2.1 

FO 

4 

17.1 

16.8 

0.15 

1.7 

F5 

5 

17.1 

20.9 

0.08 

1.3 

GO 

6 

26.4 

21.6 

0.22 

1.1 

G5 

7 

23.9 

22.6 

0.35 

0.93 

KO 

8 

19.8 

19.7 

0.34 

0.78 

K5 

9 

25.0 

25.1 

0.85 

0.69 

MO 

10 

17.3 

24.7 

1.29 

0.47 

M5 

11 

23.3 

24.1 

6.39 

0.21 

white dwarf 

12 

38.3 

36.6 

0.72 

0.9 

giant 

13 

21.0 

23.7 

0.06 

4 










